A&A manuscript no. 

(will be inserted by hand later) 



Your thesaurus codes are: 
04 (02.08.1; 03.13.4; 11.06.1) 



ASTRONOMY 
AND 
ASTROPHYSICS 
1.2.2008 



An SPH code for galaxy formation problems 

Presentation of the code 

John Hultman and Daniel Kallander 

Astronomiska observatoriet i Uppsala, Box 515, S-751 20 Uppsala, Sweden 
e-mail: pohlman@astro.uu.se 

Received 22 October 1996 / Accepted 12 February 1997 (A&A 324, 534) 



Abstract. We present and test a code for two- fluid 
simulations of galaxy formation, one of the fluids be- 
ing collision-less. The hydrodynamical evolution is solved 
through the SPH method while gravitational forces are 
calculated using a tree method. The code is Lagrangian, 
and fully adaptive both in space and time. A significant 
fraction gas in simulations of hierarchical galaxy formation 
ends up in tight clumps where it is, in terms of compu- 
tational effort, very expensive to integrate the SPH equa- 
tions. Furthermore, this is a computational waste since 
these tight gas clumps are typically not gravitationally 
resolved. We solve this by merging gas particles in these 
regions, thereby limiting the SPH resolution to the gravi- 
tational force resolution, and further speeding up the sim- 
ulation through the reduction in the number of particles. 
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Smooth particle hydrodynamics (SPH) is a fully La- 
grangian numerical method for gas dynamics. SPH was in- 
troduced by Lucy (1977) and Gingold & Monaghan (1977) 
to avoid the limitations of grid-based methods. When cou- 
pled with a fast algorithm to calculate gravitational forces, 
this method has been used successfully to study hierar- 
chical galaxy formation (Evrard 1988, Hernquist & Katz 
1989, Navarro & White 1993). 

In this paper we present an SPH implementation, cou- 
pled with a gravitational tree method, that has been 
specifically developed to study the formation of galaxies. 
The code is applied to several test cases to clarify the ad- 
vantages and limitations of the implementation. The code 
used for the simulations described in this paper is fully 
Lagrangian, does not put restrictions on the geometry of 
the problem, and has a locally varying resolution both in 
space and time. 



1. Introduction 

Recently, there has been a dramatic increase in extra- 
galactic observational data through, e.g., the Hubble and 
COBE satellites, and very large ground based telescopes. 
These detailed observations have put new demands on the- 
oretical models of galaxy formation. 

Gas dynamics with proper energy dissipation is of fun- 
damental importance in galaxy formation. Analytic mod- 
els of gas dynamics are typically restricted to systems 
possessing a high degree of symmetry. Hierarchical galaxy 
formation, on the other hand, is a highly inhomogencous 
process. In order to follow the three-dimensional evolution 
of the baryonic component in a proto-galaxy, numerical 
techniques are required. 

Send offprint requests to: J. Hultman 



Simulating the formation of galaxies in a hierarchi- 
cal model is a demanding task for any numerical method. 
Small objects form first, and the characteristic mass of 
objects then grows rapidly with time, as smaller objects 
merge into larger ones. At any given time there is a wide 
range of object masses, and the code must be able to han- 
dle many different scales simultaneously. When galaxies 
merge, part of the gas falls to the center of the new galaxy. 
This gas concentration effect is especially pronounced in 
simulations without star formation, where the baryonic 
mass remains gaseous throughout the simulation. The gas 
cores that collect at the center of dark matter halos are 
typically to small to be gravitationally resolved, but can 
still be the cause of much of the calculational cost. We 
avoid this problem by merging particles in unresolved gas 
cores, thereby limiting the hydrodynamical resolution to 
the gravitational resolution, and furthermore speeding up 
the calculations significantly. 
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2. The SPH method 



2.2. Hydrodynamical equations 



2.1. SPH principles 

The central concept in SPH is the introduction of the 

interpolation procedure. This is a procedure to define a 
continuous, differentiable function, given values at dis- 
crete points in space. We shall only outline the procedure, 
since the method has been presented many times else- 
where (Hernquist & Katz 1989, Benz 1990 and Monaghan 
1992). 

A common way to smooth a continuous field is to con- 
volve it with a suitable kernel W{x] h), 



{f{r))= j f{r')W{r-r';h)d^r', 



(1) 



the integration being over all space. The kernel, or win- 
dow function, W {x. h) goes to zero for large values of x/h. 
the parameter h thus specifies the width of the smooth- 
ing window. Furthermore, to preserve the overall normal- 
ization of the field that is smoothed, and to ensure that 
lim?i^o (/('")) = the kernel must satisfy 



/ 



W{r;h)d^r = l. 



(2) 



In SPH, particles are assigned local values of thermo- 
dynamic quantities, and corresponding continuous fields 
are defined by convolution of particle properties with a 
suitable kernel. The volume integration in Eq. (1) then 
turns into a sum over all particles 



N 



{fir)) = ^f{rj)Wi\r-r,\;hj)-^, 
where N is the total number of particles, r 



(3) 



3 ' 



ij, hj, pj 



are the position, mass, smoothing length, and density, re- 
spectively associated with particle j. The factor mj/ pj is 
the volume element that SPH associates with particle j. 

Much of the usefulness of SPH derives from the ease 
with which gradients of fields can be calculated. Ignoring 
surface terms, a direct differentiation of Eq. (3) yields 



N 



{'^f{r)) = Y.f{rj)VW{\r-rj\;hj)-l. 



(4) 



We see that once the gradient of the chosen kernel has been 
calculated, calculations of field gradients proceed similarly 
to the calculation of field values. It is especially worth 
noting, that the sets of particles that contribute to the 
SPH sums in Eq. (3) and (4) are identical. 

As is common, we use the spline kernel, that was first 
proposed by Monaghan & Lattanzio (1985), 



Kiw-l)w' + l 



< M < 1 

1 < u < 2 
u>2 



(5) 



where u = r/h and ct is I/tt in the three dimensional case. 



The time evolution of the fluid is determined by the usual 

set of hydrodynamic equations. Expressed in Lagrangian 
form, the SPH expression for the momentum equation has 
the form 



"dF 



N 



Pz P, 

4 + ^ + n. 



(6) 



where v is the velocity field and $ is the gravitational po- 
tential. Wij is the kernel used for the interaction between 

particle i and j, and H^ is an artificial viscosity pressure 
term, needed to handle shocks in SPH, described in detail 
later. For work with galaxy formation, the gas is usually 
assumed to be ideal. The pressure can then be expressed 
in terms of the density, and the specific internal energy u, 



P = (7 - l)fm, 



(7) 



where 7 is the specific heat ratio, being 5/3 for a mono- 
atomic gas, and the density is obtained through 



N 



Pi = ^mjW{\ri 



(8) 



The energy equation can be written as 



N 

E 



+ 



H,: 



r- A 



(9) 



where F and A are non-adiabatic heating and cooling 
terms, respectively. 

In order to get fully symmetric inter-particle forces, 
the kernel term Wij, must be symmetric. We have chosen 
to symmetrize in the smoothing lengths, so that Wij = 
W{rij,hij), where nj = \ri - rj\ and hij = (hi + hj)/2. 

In SPH shocks are handled by introducing artificial 
viscosity to smear out discontinuities to resolvable scales. 
We use the artificial viscosity introduced by Monaghan 
& Gingold (1983) and Monaghan & Varnas (1988). The 
viscosity term in Eq. (6) is given by 



Uij = — {-aiiijCij + (3n^ij) 

Pij 



(10) 



where pij = {pi+pj)/2, and = (cj + Cj)/2 is the average 
speed of sound of particles i and j. jXij is given by 



P-ij 



t1 ' t1 ' 



< 



ry > 



(11) 



where Vij = Vi — Vj, hij = (hi + hj)/2, and rj is a constant 
that prevents singularities for small particle separations. 
Typically, 77 = 0.1, and a and /? are set to values close to 
unity. 
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There is a disadvantage with this type of artificial vis- 
cosity, in that it introduces shear viscosity into the flow, 
in addition to the (desired) bulk viscosity effects. We cor- 
rect for this as shown in Navarro & Steinmetz (1996). All 
simulations in this paper uses this correction, except the 
Collapse test, Sect. 4.2, where essentially no shear flow is 
present. 

The smoothing length for a particle i is set by requiring 
that the number of particle neighbors, within a distance 
of 2hi, are roughly equal to a constant number, Ug. The 
smoothing length of a particle is adjusted at each time- 
step according to this criterion. To fulfill this, the smooth- 
ing length is updated at the beginning of each time-step 
according to 



1 



1/3 



(12) 



where /i* 



hi, and 



are the previous and present 



values of the smoothing length and the previous num- 
ber of neighbors respectively. The number of neighbors is 
then calculated using the predicted value of the smooth- 
ing length, and if found to deviate from Ug by more than 
ritoi, the above formula is reiterated. We chose rig = 64 
and ntoi = 10. 

We do not (yet) include the, so called, V/i terms, 
which have been closely examined by Nelson & Papaloizou 
(1993) and Serna et al. (1995), due to the relative com- 
plexity of these terms. 

2.3. Gravitational forces and neighbor finding 

We have chosen the hierarchical tree approach when cal- 
culating gravitational accelerations. It is a common choice 
together with SPH, because, like SPH, it uses no grid and 
shares a similar kind of flexibility. Moreover, the gravita- 
tional force calculation can then be efficiently combined 
with the SPH calculations. 

We use the Barnes-Hut algorithm (Barnes & Hut 1986, 
Hernquist 1987, and Hernquist & Katz 1989), which uses 
an octal tree to represent a hierarchical subdivision of 
space into cubes. We have introduced a small modifica- 
tion in the tree construction, in that cells are recursively 
subdivided until the sub-cells contain eight or less parti- 
cles. (The standard Barnes-Hut algorithm continues sub- 
dividing space until a cube contains one or zero particles.) 
This modification significantly reduces the number of cells 
in the tree, and has the additional advantage that particles 
can be placed on top of each other. The latter is sometimes 
useful when setting up initial conditions for simulations 
that involve both a gas and a dark matter component. 

When computing the acceleration or potential of a par- 
ticle, we use the standard Barnes-Hut criteria. The tree is 
traversed, root cell first, level by level, and a cell is ac- 
cepted as a force term if s/d < 9, where s is the size of 
the cell and d is the distance to the center of mass of the 
cell. 9 is usually in the range 0.5 — 1.2. Rejected cells are 



sub-divided further. We avoid self-forces by rejecting ceUs 

that the particle itself resides in. (Salmon & Warren 1993). 
Quadrupole moments are optionally included. 

The tree construction, as well as the traversal is fully 
vectorized. The "breath-first" approach was preferred over 
the "depth-first" method, because it leads to better vec- 
torization with individual particle time-steps (Makino 
1990). 

The gravitational softening is implemented with ker- 
nel softening (Gingold & Monaghan 1977 and Hernquist & 
Katz 1989); where the force between two particles is calcu- 
lated by considering one particle to have a radial density 
distribution given by the same spline kernel as for the SPH 
interactions, Eq. (5). 

Gravitational softening parameters ej can be specified 
for each particle. Typically, the softening will be different 
for the dark matter component and the gas component, 
and sometime regions with different resolution are used. 
The softening is chosen so as to keep two-body relaxation 
at reasonable levels. Typically, we calculate the two-body 
relaxation time-scale (Binney & Tremaine 1987) for grav- 
itationally bound objects as 



Pretax — 



NR 



8{v) log A' 



(13) 



where N, R and (v) are, respectively, the total num- 
ber of particles, radius and mean velocity, for the object 
under consideration, log A — log{bmax / bmin) , is known 
as the Coulomb logarithm. We use impact parameters, 
bmax = RN'^^^, according to Smith (1992). bmin = e-min, 
the minimum softening parameter in the object. 

The need to symmetrize pairwise forces occurs also 
here, and for particle-particle interactions we symmetrize 
the softening lengths, = (e^ + ej)/2, in the same way as 
the SPH smoothing lengths. The particle-cell interactions 
are not straight-forward to symmetrize in a hierarchical 
tree method, and we do not do this. 

Instead of starting from scratch, we begin the neighbor 
search where the tree traversal for the gravitational force 
stops, making the calculations more efficient. The gravity- 
neighbor calculations typically constitutes over 60% of the 
computing time. The gravity-neighbor package is calculat- 
ing for one particle at a time, and can therefore be exe- 
cuted in parallel. It was straight forward to implement this 
on shared memory (SMP) architectures, like the CRAY Y- 
MP. Generally this reduces calculation times by a factor 
of two. 

2.4- Time integration 

We use the second order Runge-Kutta (RK) integra- 
tor, with individual time-steps, described by Navarro and 
White (Navarro & White 1993). The main reason for this 
choice is that it allows the time integration errors to be 
estimated, and thereby gives good control over particle 
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time-steps. There is no need to directly incorporate any 
stability criteria, e.g., the Courant condition. 

The system is always advanced with the smallest time- 
step, and a subroutine keeps track of which time- levels of 
particles are at either the beginning of a time-step, or at 
the middle of a time-step, and need to have derivatives 
of dynamical variables calculated. A linked list structure 
is used for the bookkeeping of the particles different time 
levels. The overhead caused by this book-keeping consti- 
tutes less than 1% of the computing time for a typical 
simulation. 

We use an implicit scheme for the integration of the 
energy equation, similar to Hernquist & Katz (1989) and 
Anzer et al. (1987). A trapezoidal formula 

u^+i = + (u^+^ + OAii/2, (14) 

is solved iteratively, together with Eq. (9), at each indi- 
vidual time-step At^. A combined Newton- Raphson and 
bisection scheme, as described in Press et al. (1992), is 
used to find the root of the energy equation (14). In ab- 
sence of cooling and heating terms, Eq. (9) is linear in Uj, 
and Eq. (14) is solved in one iteration cycle. 

When strong radiative cooling is present, the integra- 
tion tolerance criterion on the internal energy equation 
would become very prohibitive, regardless of what type of 
criterion is used, (Courant e.t.c). Physically, this occurs 
when the cooling time-scale becomes much shorter than 
the local dynamical time-scale 

It would be extremely expensive in terms of computing 
time to integrate the system at the time-step imposed by 
the cooling time-scale. Moreover, for the purpose of track- 
ing the dynamical evolution of the system, it is not nec- 
essary to let the time-step become much shorter than the 
local dynamical time-scale. Therefore, the accuracy con- 
straint on the energy equation is not allowed to limit the 
time-step below 10~^ • td.yn- 

Although the implicit scheme is unconditionally sta- 
ble, there will still be errors in the time integration. We 
have found that a damping of the cooling terms, simi- 
lar to Katz & Gimn (1991), improves the accuracy of the 
time integration, without affecting the dynamical evolu- 
tion. A particle is not allowed to lose more than half its 
internal energy during one time-step, by cooling. We use 
a sharp cutoff on the cooling part of the derivative. This 
has the effect of broadening the fast time variation into a 
resolved exponential decay over a few time-steps. The rea- 
son why this works is that it occurs only in regions where 
the cooling time-scale is much shorter than the dynam- 
ical time-scale. The energy that is not radiated away in 
one time-step, will be radiated away in the following few, 
while the other variables, like density and velocity, can be 
considered to be roughly constant. 



Typical radiative cooling functions often have multiple 

peaks. We have noted this may cause miiltiple roots in the 
solution of the energy equation. Of course, only one root 
is correct, i.e., the one closest to the root at the previ- 
ous time-step. We solve this problem by finding all roots 
and accepting only the closest one. This may seem both 
computationally costly (and mathematically impossible, 
but the cooling functions used, are tabulated at a finite 
number of points. Typically, 128-256, and therefore these 
calculations take insignificant computational time. 

3. Merging 

3.1. Merging of particles 

Our code was written with simulations of hierarchical 
galaxy formation in mind. A typical feature of such sim- 
ulations is extreme clumping of matter. The gas collapses 
into very dense gravitationally bound objects, where cool- 
ing time scales are much less than the local dynamical 
time-scale. Some of these objects contain very little angu- 
lar momentum, and the gas therefore continues to contract 
until the contraction is artificially halted by the finite nu- 
merical force resolution. This clumping leads to two prob- 
lems. 

First, the suppression of gas contraction on unresolved 
scales is completely artificial, and also puts a limit on the 
gravitational binding energy. This makes it easier for gas 
to be ejected from bound objects during violent mergers, 
and could lead to an underestimate of the fraction of mat- 
ter that has collapsed to galactic densities at a given time. 

Second, this is a challenge for SPH because accurate 
time integration of the gas in extreme high density re- 
gions requires large amounts of computational work. We 
have chosen to set the SPH smoothing lengths by requiring 
a roughly constant number of interacting neighbors. This 
leads to very small smoothing lengths in regions with ex- 
tremely high density. To maintain stability, it is therefore 
necessary to use very small time-steps in these regions. We 
have found that the computational cost for evolving gas 
in gravitationally unresolved cores of compact objects can 
dominate the total run-time. This is a very undesirable 
situation. 

The problems with small time-steps in unresolved re- 
gions can be somewhat alleviated, by requiring that the 
SPH smoothing always be larger than the gravitational 
smoothing. This limits the maximum gas resolution to 
roughly the gravitational force resolution. However, the 
problem then occurs in the SPH force evaluations instead. 
Setting a lower limit on the SPH smoothing in high density 
regions leads to large numbers of interacting neighbors. 
This makes the force evaluations very expensive, due both 
to the large amounts of terms in the SPH loops and the 
cost of finding the SPH neighbors. The code deteriorates 
towards an inefficient implementation of direct summa- 
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tion. Limiting the hydrodynamical resolution in this way 
actually slows the code down. 

The reason that the run-time can increase when the 
resolution is decreased, by setting a lower limit on the SPH 
smoothing length, is that the accuracy of the SPH force 
calculations increases. A lower limit on the SPH smooth- 
ing length leads to large numbers of hydrodynamically in- 
teracting neighbors in unresolved regions, and therefore 
reduced sampling errors in the SPH force calculations. 
An alternative would be to only use a random subset of 
the interacting neighbors in the SPH calculations. This 
is consistent with the SPH formalism, but it introduces 
an extra source of noise in forces. It would also be inef- 
fective to first find all neighbors of a particle, and then 
discard most of them only to use a small subset in the ac- 
tual force evaluation. A tree search can be implemented to 
find a pseudo-random subset of a particles neighbors, one 
by one, by using a depth-first search. This would make it 
unnecessary to find large numbers of neighbors. It looks 
non-trivial to vectorize such a scheme, and we have found 
that large subsets have to be used for good accuracy. 

We avoid all these problems by not imposing any lower 
limit on h, thus ensuring efficient force evaluations, and 
by instead merging particles in high density regions, thus 
avoiding small time-steps and further speeding up the sim- 
ulation by reducing the number of particles. This has pre- 
viously been tried with success by Monaghan & Varnas 
(1988) for simulations of cloud collapse, and by Steinmetz 
(1996) for galaxy formation problems. 

Merging of particles in high density regions can keep 
the hydrodynamical resolution from exceeding the gravi- 
tational resolution, without slowing down the force evalu- 
ations, thus greatly speeding up simulations without sac- 
rificing accuracy. 

The first step is to identify high density regions 
that are severely affected by gravitational smoothing. In 
these regions it is pointless to calculate substructure on 
scales smaller than the gravitational smoothing length. We 
therefore merge particles in these regions into more mas- 
sive, and fewer, particles. The dynamics on these small, 
unresolved, scales cannot be followed in our code, with 
or without merging of particles. The dynamical effects of 
particle merging should be small outside the regions where 
merging occurs. To satisfy this requirement all global dy- 
namical quantities must be left unchanged by the merging 
process. Thus, when merging two particles into one, global 
quantities like mass, kinetic energy, internal energy, poten- 
tial energy and angular momentum must be preserved. 

To make sure that global dynamical quantities are ac- 
curately conserved, the following five conditions have to 
be met before two particles are allowed to merge. ( Ai , A2 , 
etc. are tolerance parameters. The quantities rrii, pi, Vi, 
Ci and (f)i are the mass, local gas density, velocity, grav- 
itational softening length and gravitational potential of 
particle i. A "cm" subscript indicates a quantity that is 
evaluated in the center of mass of the two merging par- 
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tides. Vij = Vi — Tj is the relative position of the two 
merging particles.) 

I) The two particles are in a region where gravitational 
forces are severely affected by limited force resolution. 

min (pj,Pj) > ^2 niax(pei,Pe,) (16) 
where 

= 32;^' ('^^ 

which is the density that would be calculated for parti- 
cle i if all particles were evenly distributed within radius 
2ei, and they all had mass m^. In regions where the den- 
sity is much higher than p^, the mean inter-particle dis- 
tance will be much less than e, and gravitational forces 
will be affected by the limited resolution. Using ngmi, in- 
stead of summing up the mass of all neighbors, makes 
the density criterion stricter for more massive particles, in 
regions with particles of varying masses. This evens out 
spatial particle mass fluctuations, by letting less massive 
particles merge before more massive ones. 

II) The two particles are SPH neighbors and within 
one smoothing length of each other. 

The particles should be close together to minimize the 
artificial displacement of gas properties that occurs when 
two particles are merged. 

HI) The resultant particle mass of the merger is below 
the preset limit 

rui < Airriinitiai, (18) 

where rriinitiai is the gas particle mass before any merging 
has taken place. This is a limit on the amount of merging 
that can occur in a simulation, and can be useful as a 
guarantee that the gas resolution does not fall below a 
preset limit. 

IV) The angular momentum lost is small. 

I'^ij X (Pic„ -Pjcm)\ < A^Lscale (19) 

where L scale is a relevant angular momentum scale for the 
problem at hand, which in all cases discussed here is equal 
to the total angular momentum of the system. 

V) The kinetic energy lost is small compared to the 
gravitational potential energy, 

"^^^ {vi - Vjf < Aiirrii + mj)4)i. (20) 

VI) The particles are synchronized in time and have 
just completed a full time-step. This is necessary to avoid 
breaking the time integration scheme. 

These criteria are ordered so that the least compu- 
tationally expensive ones are checked first. The compu- 
tational cost for the whole merging procedure is then 
negligible. Whenever criteria I- VI are met the two par- 
ticles merge, and form a new particle. The new particle is 
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given properties to conserve mass, momentum and ther- 
mal energy. Furthermore, we let the gravitational smooth- 
ing be proportional to (m^ + nij ) 3 . Changing gravitational 
smoothing changes the potential. If conservation of bind- 
ing energy is desired, it is safer to leave the gravitational 
smoothing unchanged. The reason is that the internal en- 
ergy have to be changed and this may cause unphysical 
decrease in entropy. 

When two particles merge, the resolution in that re- 
gion changes. This is accompanied by a small change in the 
SPH smoothing length in the region. To avoid large dis- 
continuous changes in SPH quantities, which could lead to 
spurious effects in the time integration scheme and could 
set up artificial shock waves, a given mass clement is not 
allowed to merge more than once per local dynamical time. 
This ensures that particle merging is a smooth process 
where SPH quantities have time to adjust to changes in 
the resolution. We have found this criterion to be crucial 
to accurately reproduce the results of simulations without 
particle merging. 

Throughout this paper the values of the tolerance pa- 
rameters are Ai = 10,^2 = A3 = A4 = 10~*. With 
these vahies we have always found that energy and angular 
momentum are conserved to well within one per cent. It 
should be noted, this choice of parameters allow for rather 
aggressive merging, sometimes making more than half the 
gas particles to merge, but all merging occurs in unre- 
solved cores, and it does work well in the tests presented 
here. 



4. Tests 

4.1. Test of the gravitational interaction; King models. 

For objects like galaxies, which have a tremendous num- 
ber of constituent bodies, the dynamics may be described 
by the Boltzmann equation. A dark matter component 
can usually be considered collision-less in contrast to the 
collisional gas component, for which the Boltzmann equa- 
tion can be well approximated with the hydro dynamical 
equations. In an N-body method discreetness effects are 
unavoidable, due to the finite number of particles used. 
A practical way to test an N-body code, is to run it on 
a steady state system. Conservation of energy, and lin- 
ear and angular momentum, can be examined. Also im- 
portant is that the effects of two-body relaxation may 
be investigated. For a steady state system the potential 
is time-independent, and thus the energies of individual 
particles should be conserved. Moreover, if the system 
is spherically symmetric, individual particle angular mo- 
menta are conserved. Due to the motions of the particles 
in an N-body model the potential will never be strictly 
time-independent, the forces not strictly central, and there 
will always be some relaxation. The particles energies and 
angular momenta will diffuse, and this effect can be mea- 
sured. 



The most practical static models to choose for testing 

purposes are King models (King 1966), since they are of 
finite extent in phase space. Following Binney & Tremaine 
(1987), King models have a phase space density 



MS) 



pi(27r(j2)-i(e^/<^'-l) 




£>0; 
£ <0, 



(21) 



where the relative energy, £, and the relative potential, ^, 
are defined as 



£ = -E + ^0 = - -v^ and * = -h $0 



(22) 



The parameter a is a measure of the one dimensional 
velocity dispersion, and pi is a mass-normalization con- 
stant. Integrating over velocity space we get the density 
as a function of ^, 



PKi'S) = Pi 




(23) 



where erf{x) is the error function. ^ will satisfy the Pois- 
son equation, 



dr 



dr 



(24) 



which is an ordinary differential equation for 5'(r) and can 
be integrated numerically once we have suitable boundary 
conditions. For positive £, which is our region of interest, 
we must have positive , thus the value of Vl/(0) is one 
condition. We note that the enclosed mass inside radius r 
can be given in terms of as 



M(r) = - 



r^d* 
G d^' 



(25) 



Restricting ourselves to a finite central density, we see 
that {d'^ /dr) = at r = is a natural condition. Since 
^"(0) is positive and {d^ /dr) is negative for r > 0, \E'(r) 
must decrease and become zero for some radius rt where 
also the density will vanish. This radius is known as the 
tidal radius. 

King models are commonly parameterized in terms of 

^'(0)/cr^, because the dispersion parameter a just spec- 
ifies the velocity scale. The mass scale is still free, but 
can be specified with a change of scale, i.e., by adjusting 
the parameter pi. The radial Poisson equation is solved 
numerically using a Runge Kutta method. The enclosed 
mass M(r), is calculated in order to distribute the parti- 
cles. Particles arc initially placed on a grid in a spherical 
region and then the grid is stretched into the desired den- 
sity profile. 

In general, it can be hazardous to use grid-setups for 
very cool collision-less systems, but this will not be the 
case here since the particles are started with substantial 
velocities that are randomly distributed. The advantage 
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Fig. 1. Density profile of a King model with central potential 
^((0) = 5(7^. Units are given by G = M = i? = 1. 



with the grid placement, over a random placement, is that 
random fluctuations are minimized, and it is easier to com- 
pare the realization with the analytic potential. The sys- 
tem will also be closer to equilibrium. 

Particles arc given random velocity directions, isotrop- 
ically distributed. Random velocity magnitudes are then 
assigned to the particles, with a distribution given by Eq. 
(21). 

The velocity probability density is given by 



fv{v) = fK{v\r) = 



(26) 



where fK{v\r) is the conditional probability density and 
Mr) = J fK{r,v)d\ = pK{r). (27) 



The distribution function is then 



V 

Fv{v)=A'K j JkW 





r)v'^dv', 



(28) 



and can may integrated in closed form. 

In order to be able to compare results with the tests of 
Hernquist & Barnes (1990) and Huang et al. (1993) we set 
up a King model with central potential ^'(0) = Scr^. Units 
were chosen such that the gravitational constant G = 1, 
the total mass M = 1 and dispersion parameter a = 0.762. 
This gives a total three-dimensional velocity dispersion of 
unity, and thus a total energy of —1/2. The tidal radius is 
2.18. 

Four different runs were made with N = 4626 or 
N = 15238. A softening parameter e = 0.025 was used 
in all runs, and they were continued to time t = 8, which 
corresponds to around 20 half-mass dynamical times. The 
Barnes-Hut parameter 6 was 1.0 in all runs, except run 



3, where ^ = 0.1 was used. All runs used individual time- 
steps, except run 4, where a constant time step At = 1/48, 
was used. In order to test the relaxation effects we exam- 
ined the relative changes in particle energies 



x, = {E,{h)-E,{t2))/E,{t2), 



(29) 



between times t\ and where Ei = v1/2 + $j In Fig. 2, 
the results for run 1 is shown. In Table 1, the standard 
deviation S and mean value x are presented for the runs. 
Absolute deviations were also examined, but they had the 
same behavior as the standard deviations. The changes 
in energies are expected to be a random walk diffusion 
process because of small random accelerations, caused by 
the random noise in the particle forces. The diffusion is 
expected to behave such that 



(30) 



where Dg is an empirical diffusion constant. Ds was cal- 
culated with a linear least squares fit of log 5 as function 
of log At. The diffusion rates are slightly larger, but com- 
parable to those of Hernquist & Barnes (1990). 

Conservation of energy, linear and angular momentum, 
and center of mass position was also examined. The max- 
imum deviations during each run, are listed in Table 1. 
The time evolution of 5 is plotted in Fig. 3. The straight 
line implies a fairly good fit of Eq. (30) . In agreement with 
what was found by Hernquist & Barnes (1990), the only 
parameter having a significant effect on the diffusion rate 
was the number of particles used. 

Judging from the CPU time used, run 4 seems to be 
over two times as fast as run 1, but run 4 uses a constant 
time step, five times larger than the smallest time step 
of run 1 . Run 1 conserves energy roughly a factor of four 
better than run 4. Three time levels are occupied for run 
1, and the benefit of individual time-steps is roughly a 
factor of two. 

4-2. The collapse test 

One of the most common tests for astrophysical SPH 

codes is the collapse test of Evrard (1988). It has also 
been presented by Hernquist & Katz (1989), Steinmetz & 
Miiller (1993), Nelson & Papaloizou (1994), and Serna et 
al. (1995). One dimensional, (spherically symmetric), fi- 
nite difference solutions have been calculated by Thomas 
(1987) and Steinmetz & Miiller (1993). We will henceforth 
refer to this test as the "collapse test". The reason why 
this test is popular is that it is simple to set up, and that it 
tests an "SPH + gravity" code on the aspects of adiabatic 
flow, shocks and gravitational collapse. The initial setup 
for this problem is a gas sphere, of mass M and radius i?, 
with density profile 



p{r) 



M 1 



27ri?2' 



(31) 



8 John Hultman et al.: An SPH code for galaxy formation problems 

Table 1. Data for the various King-model runs. Relative changes in particle energies, between times t = 4 and t = 8. Conserved 
quantities. 



Run 


N 


e 


At 


s 


X 


Ds 


AE/E 


AL/L 


ARcM/n 


AVcm/V 


CPU 


1 


4626 


1.0 


0.0039* 


0.078 


-3.3 X 10"^ 


0.035 


5.4 X 10"^ 


4.8 X 10"" 


0.011 


0.036 


1.94 


2 


15238 


1.0 


0.0063* 


0.041 


-7.8 X 10"" 


0.020 


5.5 X 10"'^ 


1.4 X 10"'' 


8.1 X lO"-'* 


0.026 


9.39 


3 


4626 


0.1 


0.0039* 


0.080 


-5.4 X 10"^ 


0.036 


5.6 X 10"^ 


0.016 


2.3 X 10"^ 


7.3 X 10"'^ 


32.1 


4 


4626 


1.0 


0.020 


0.076 


-0.011 


0.036 


0.021 


7.2 X 10"" 


5.8 X 10"^ 


0.019 


0.83 



* shortest individual time step. Scales L and V are of order unity. CPU- times are in hours on a HP 735 workstation. 



-3.5 -3.0 -2.5 -2.0 -1.5 -1.0 
Total energy El (Internol units) 



-0.5 0.0 



Fig. 2. The relative change in individual particle energies, as 
a function of the energies, Measured at times t = 4 and t = 8. 
Internal units are given byG = M = i? = l. 
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Fig. 3. The standard deviation of relative particle energies, S, 
as a function of the time, for run 2. The slope of the dashed 
line is 0.506. 



The gas is initially isothermal with an internal en- 
ergy u = 0.05GM/R, and the velocity is zero everywhere. 
Other test parameters are the specific heat ratio 7 — 5/3, 
the artificial viscosity parameters a = 1, P = 2 and 
r] = 0.1, the Barnes-Hut tolerance for the gravitational in- 
teraction ^ = 0.8 and a softening parameter e = 0.1. The 
number of particles used was N = 3828, and the number 
of neighbors in the SPH summations was =64. 

We use the same kind of stretched grid set up, as in 
Sect. 4.1. It gives a more relaxed initial configuration than 
the corresponding random distribution. For this particular 
distribution the errors in the SPH-density are very small, 
as can be seen in Fig. 5. 



Time: 0.0000 (Internal units) 



Rodius (Internol units) 



Fig. 5. The initial l/r-density profile for the collapse test. 

The dotted line is a plot of the desired density, Eq. (31). The 
SPH density of the particles follows the desired density, except 
at the very center r < 0.04, and at the outer edge r > 0.9, 
where boundary effects are visible. Internal units are given by 
G = M = R = l. 



A good way to check for SPH sampling errors, is to per- 
form an SPH evaluation of a constant unity field, putting 
f{rj) = 1 in Eq. (3). It is important that no edges etc., 
where the sampling necessarily is bad due to the neglect 
of surface terms are within the investigated volume. Cor- 
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Fig. 4. Density, pressure, internal energy, velocity and Mach number, for some different times, in the "collapse test". Units are 
given hyG = M = R=l 
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respondingly, the gradient of the unit field may be exam- 
ined. The statistics of the SPH representations of a unit 
field and its gradient, for the initial configuration in Fig. 
5, are presented in Table 2. Although the result for a unit 
field shows very little spread for this particular configura- 
tion, the result for the corresponding gradient field shows 
considerably more spread. However, these fluctuations are 
much smaller than those corresponding to a random place- 
ment setup, as shown in Table 3. To avoid boundary ef- 
fects, when calculating the numbers in Table 5 and Table 
3, only a selection of 1800 particles within radius 0.7 were 
used. 



Table 2. SPH sampling errors for a 1/r-density profile, using 
a stretched grid setup. 



quantity 


mean 


max 


min 


a 


unit field 


1.000 


1.002 


0.9989 


4.495 X 10"" 


|V(uuitfi(-ld) 


0.07012 


O.i!)."):! 


O.OliOT 


().0:',()0i 



Table 3. SPH sampling errors for a 1/r-density profile, using 
random placement. 



quantity 


mean 


max 


min 


a 


unit field 


1.008 


1.352 


0.6744 


0.1024 


|V(unitfield)| 


2.898 


9.989 


0.2110 


1.3079 



Units in the collapse test are chosen so that G = M = 
R = 1, and the usual plots of density, pressure, internal 
energy, velocity and Mach number are shown in Fig. 4. 
The results show good agreement with what have been 
reported by the authors mentioned above. 



Table 4. Conserved quantities for the collapse test. 



TV AE/E 


AL/L 


ARcm/R 


CPU 


3828 9.7 X 10"^ 


1.0 X 10"* 


3.8 X 10-* 


5.25 



Scale L are of order unity. CPU-time are in hours on an HP 
735 workstation. 

With only 4000 particles the outgoing shock is not 
very sharp, but as noted in Steinmetz & Miiller (1993) 
all global quantities are still reproduced rather well. Usu- 
ally the smearing of shock fronts is rather severe in SPH, 



2 - 




-3 L : : : : : : : : I : : : : : : : : : I : : : : : : : : ^ 

12 3 

Time 

Fig. 6. Energies plotted against time for the collapse test. 
From top to bottom are the curves for the thermal, kinetic, 
total and potential energy. Units are given hy G = M = R = 1. 

but still this affects the evolution of the shock fronts sur- 
prisingly little. Depending on the problem, one often needs 
some factors of 10000 particles in shock regions, to see the 
shocks clearly resolved. As for the King-model test, some 
conservation properties were examined and are listed in 
Table 4. The shortest (individual) time-step during the 
simulation was 10"**, and the particles occupied a span of 
7 time-levels. In Fig. 6, the total thermal energy, total ki- 
netic energy, total potential energy and total energy, are 
plotted. The maximum thermal energy is slightly lower, 
than reported by Hernquist & Katz (1989) and Steinmetz 
& Miiller (1993). This is because of the higher number 
of neighbors used, here being 64. as compared to the 32 
neighbors used by above mentioned authors. 

4-3. Test with strong radiative cooling 

As a test including radiative cooling, and merging of par- 
ticles, we have simulated the collapse of a rotating sphere, 
with parameters reminiscent of a proto-galaxy. This test 
case has been studied by Navarro & White (1993) and 
Serna et al. (1995). 

The initial density field is spherically symmetric, with 
a radial density profile p oc 1/r. The sphere is started in 
solid body rotation with a dimensionless spin parameter 
A = J|£;|i/2/G'Mf/i^ « 0.1 (J and E are the total angu- 
lar momentum and total energy). The initial radius is 100 
kpc, and the total mass is lO^^M©, with 10% of the mass 
in gas and 90% of the mass in dark matter. The gas starts 
at a temperature of 1000 K. The gas and the dark matter 
components are represented by 2000 particles each. Grav- 
itational softening parameters were taken to be 2 and 5 
kpc for the gas and dark matter, respectively. 

In these simulations a significant fraction of the gas 
falls into unresolved clumps in the central parts of the disk. 
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It is therefore a demanding test of the particle merging 
scheme. 

Fig. 7 shows the evolution of the gas and dark matter 
distributions for a simulation without particle merging. 
The dark matter virializes soon after the main collapse, 
whereas the gas forms a thin disk. Collisional line cool- 
ing is almost completely effective in radiating away ther- 
mal energy from adiabatic compression and shock heating. 
Only a slight fraction of the gas is heated to temperatures 
above 30,000 K (see Fig. 10), and thermal pressure does 
not play any significant role in the evolution. Radial ve- 
locities are effectively dissipated, and most of the gas has 
formed a rotationally supported disk after one collapse 
time. 

The disk displays clear spiral structures after t « 160. 
As reported by Navarro & White (1993), the disk is unsta- 
ble and starts breaking up into small clumps towards the 
end of the simulation. Navarro & White found that if the 
gas mass fraction is reduced to 2%, thereby increasing the 
Toomrc (1964) stability parameter, the disk is substan- 
tially stabilized. Serna et al. (1995) found the stability 
of the disk, when the gas mass fraction was 10%, to be 
intermediate to the 10% and 2% gas mass fraction case 
of Navarro & White. Our results resemble those of Serna 
et al. The reason that we get a more stable disk than 
Navarro & White, is probably due to the fact that we set 
the SPH smoothing so as to acquire w 64 particle neigh- 
bors (rig = 64), whereas Navarro & White uses roughly 32 
neighbors. When we reduce the number of SPH neighbors 
to ~ 32, we find that the disk evolution closely matches 
that of Navarro & White. 

In the simulation that includes our prescription for 
merging particles in high density regions (see Fig. 8), the 
number of particles decrease as more particles collapse 
into high density regions and merge with other particles. 
As can be seen from Fig. 9, the number of particles has 
been reduced to half the initial value at the end of the 
simulation. 




£1 

E 

500 



0.2 0.4 0.6 0.8 1 1.2 1.4 

Gyrs 

Fig. 9. Number of gas particles, in the simulation with merg- 
ing of particles, as a function of time. The number decreases 
as particles in high density regions merge into more massive 
particles. 



The gas cooling rate depends on the square of the local 
gas density. When particles are merged the local resolution 
decreases, and small scale fluctuations in the density field 
are damped. This could potentially have significant effects 
on the gas cooling. This is a fundamental problem in all 
gas dynamical simulations of galaxy formation. The gas 
density field has to be assumed to be reasonably smooth 
on scales smaller than can be resolved in the simulation. 
Fig. 10 shows that the fraction of hot gas is not altered 
by the inclusion of particle merging. 




0.2 0.4 0.6 0.8 1 1.2 1.4 



Gyrs 

Fig. 10. The mass fraction of gas with a temperature above 
3 • 10'* K for the simulation, without merging (solid line), and 
with merging of particles (dashed line) . 

The rotation curves, {vc = {GM{r)/r)i), for the sim- 
ulations are shown in fig 11. The mass distribution when 
merging is included is in excellent agreement with the 
more conventional run, without merging. The rotation 
curves are approximately flat out to the edge of the disk at 
w 25kpc. The disk that forms in the run without particle 
merging is slightly more concentrated, with a 30% higher 
gas mass inside 3 kpc. Fig. 7 and 8 indicates that particle 
merging stabilizes the disk, and slightly suppresses the for- 
mation of substructure. This is a natural consequence of 
the decreasing gas mass resolution and gravitational force 
resolution, when particles merge. 



350 




5 10 15 20 25 30 35 40 



Fig. 11. Circular velocity curves, with merging of particles 
(solid line), and without merging of particles (dashed line). 
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Fig. 7. Collapse of a rotating sphere consisting of gas and dark matter, and with a 1/r-density profile. All numbers are given 
in the same units as Navarro & White (1993). (Distances are in kpc, and the time unit is 4.71 Myrs.) 
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Fig. 8. Collapse of a rotating sphere consisting of gas and dark matter, and with a 1/r-density profile. All numbers are given 
in the same units as Navarro & White (1993). (Distances are in kpc, and the time unit is 4.71 Myrs.) Particles were allowed 
to merge during the simulation. Note that individual gas particle masses vary, and therefore that the number density of gas 
particles is not proportional to the gas mass density. 
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The total CPU time for the simulation is halved when 
particle merging is allowed, At the end of the simulation 
the CPU cost per unit time has been reduced by a factor of 
five, and the total CPU time for the simulation is halved, 
when particle merging is allowed. 

4-4- Formation of a galactic object 

In order to test our code, and especially the particle merg- 
ing scheme, on problems that arc typical of those we wish 
to study, we simulate the collapse of a proto-galaxy that 
has been set up consistently with the CDM cosmological 
model, in a fashion similar to Katz & Gunn (1991). The 
same starting conditions were used for three different sim- 
ulations, two with particle merging and one without. 



lapsed object as a function of redshift. Collapsed gas ob- 
jects were identified using a fricnds-of- friends algorithm, 
grouping particles together that had an over-density ex- 
ceeding 1000. The magnitude of the mass build-up over 
time is very similar in all three simulations. 

The mass fraction of gas with a temperature exceeding 
3 • W^K, as a function of redshift, is shown in Fig 16. The 
two 8000 particle simulations differ slightly, with roughly 
1% more hot gas being produced after z=1.2 in the simu- 
lation including merging. The curve for the high resolution 
simulation deviates clearly from the other two after z = 1. 
Between 1 < 2 < 0.5 more hot gas is produced, and after 
2; = 0.5 more hot gas is able to cool than in the lower res- 
olution simulations. At z=0 the mass fraction of hot gas 
is close to 10% for all three simulations. 



The Zeldovich approximation together with a standard 
CDM model (O = 1, A = 0, fife = 0.05, hwo = 0.5, 
as = 0.66, Bardeen ct al. 1986) power spectrum, was used 
to set up a cosmological density field. The system is given 
an initial over-density corresponding to a 3 cr peak in the 
CDM spectrum when it is convolved with a top hat fil- 
ter of mass lO^^M©. This density field was realized in- 
side a sphere with a co-moving radius of 1.46 Mpc, using 
both gas and dark matter (collision-less) particles. Gas 
and dark matter particles had a gravitational smoothing 
of 2 and 5 kpc, respectively. The system is started in solid 
body rotation, corresponding to a dimensionless spin pa- 
rameter of ~ 0.05, in an attempt to roughly approximate 
the effects of tidal interactions. The initial gas tempera- 
ture is 10,000 K, and the gas cooling rate is that of a gas 
in coUisional ionization equilibrium and with a 0.1 solar 
metallicity, as given by Sutherland & Dopita (1993). 

Three simulations were made, differing only in the 
number of particles used, and whether or not particles 
were allowed to merge. Two simulations were started with 
8000 particles, half of them used to represent the gas and 
the other half to represent the dark matter component. 
The only difference between these two simulations is that 
one employed the previously described scheme for merging 
gas particles, and the other did not. The third simulation 
was made with ten times more gas particles, and with par- 
ticle merging. The same realization of initial conditions 
were used in all three simulations, in order to make the 
final initial conditions as similar as possible. 

The systems were started at z = 30, and evolved to 
z=0. The main collapse of the proto-galaxy occurs at 
2; w 2. At z=0 a single dominant gas object with galactic 
densities has formed in all the simulations. This galactic 
object consists of a compact core surrounded by a thin 
disk, Fig 12, Fig 13, and Fig 14. This object is built up 
from a combination of continuous in-fall and the merging 
of smaller collapsed objects. 

The mass build-up of the final objects can be seen in 
Fig 15, which shows the gas mass of the most massive col- 
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z 

Fig. 15. The gas mass of the most massive collapsed object as 
a function of redshift, for 4000 particles no merging (solid line), 
4000 particles with merging (dashed line), and 40,000 particles 
with merging (dotted line). 
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Fig. 16. The mass fraction of the gas that has a temperature 
exceeding 3 • 10*/^ as a function of redshift, for 4000 particles 
no merging (solid line), 4000 particles with merging (dashed 
line), and 40,000 particles with merging (dotted line). 
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Fig. 12. Three orthogonal views of projected gas particle positions for the simulation with 4,000 gas particles and no particle 
merging. Frames sizes are 136 kpc. 




Fig. 13. Three orthogonal views of projected gas particle positions for the simulation with 4,000 initial gas particles and 
particle merging. Frames sizes are 136 kpc. Note that the individual particle mass is not constant, and that the number density 
of particles is therefore not proportional to the mass density. 



Fig. 14. Three orthogonal views of projected gas particle positions for the simulation with 40,000 initial gas particles and 
particle merging. Frames sizes are 136 kpc. Note that the individual particle mass is not constant, and that the number density 
of particles is therefore not proportional to the mass density. 
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These results seem to indicate that the effects of ap- 
plying the particle merging scheme arc small, outside the 
gravitationally unresolved cores of collapsed gas objects. 
Furthermore, a ten-fold increase of the initial number of 
gas particles in a simulation produced only moderate dif- 
ferences, lending some tentative support for low resolution 
SPH simulations of galaxy formation. 

Evolving a system with 4000 gas and 4000 dark matter 
particles initially, took 64 CPU hours on a CRAY-YMP. 
The same system with merging of particles required only 
13 CPU hours. The high resolution system required 63 
CPU hours, almost the same as the low resolution system 
without merging. The decrease in the CPU time required 
for a simulation, due to particle merging, will vary with the 
problem and the tolerance parameters used in the merging 
scheme. 
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